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Latent  heat  thermal  energy  storage  (LHTES)  has  recently  attracted  increasing  interest  related  to  the  thermal 
applications,  and  has  been  studied  by  researchers  using  theoretical  and  numerical  approaches.  The  heat 
transfer  mechanisms  during  the  charging  and  discharging  periods  of  the  phase  change  materials  (PCMs)  and 
two  basic  problems  for  phase  transformations  have  been  discussed  in  this  paper.  ID,  2D  and  3D  popular 
mathematical  solutions  based  on  the  heat  transfer  mechanisms  of  conduction  and/or  convection  have  been 
analyzed.  Then,  various  numerical  models  for  encapsulated  PCMs  in  terms  of  macro-encapsulated  PCM  and 
micro-encapsulated  PCM  were  investigated.  The  numerical  simulation  of  heat  transfer  problems  in  phase 
change  processes  is  complicated  and  the  achieved  results  are  approximate  according  to  a  number  of  conditions. 
The  accuracy  of  numerical  solution  depends  on  the  assumptions  that  are  made  by  the  authors.  Therefore  this 
review  will  enable  the  researchers  to  have  an  overall  view  of  the  mathematical  and  numerical  methods  used 
for  PCM's  phase  transformations.  It  offers  the  researcher  a  guideline  of  selecting  the  appropriate  theoretical 
solutions  according  to  their  researches'  purposes. 
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Nomenclature 

Greek  symbols 

f 

melt  fraction 

A 

difference 

c 

specific  heat  at  constant  pressure 

cc 

thermal  diffusivity  of  PCM  (m2/s) 

H 

total  volumetric  enthalpy  0/m3) 

' 1 

dimensionless  number  used  in  derivations  as  a  tern- 

h 

average  heat  transfer  coefficient  (J/kg) 

porary  substitution 

k 

thermal  conductivity  (W/mK) 

A 

dimensionless  number  in  solution  to  Neumann  pro- 

L 

latent  heat  0/kg) 

blem,  liquid  fraction 

St 

Stefan  number 

P 

density  (kg/m3) 

t 

time  (s) 

T 

temperature  (I<) 

Subscripts 

Tm 

melting  temperature  (K) 

1 

liquid 

Ts 

initial  temperature  of  solid  PCM  (K) 

s 

solid 

To 

constant  temperature  imposed  on  x=0  (K) 

i 

ith  spatial  step 

P 

present  nodal  point 

j 

jth  time  step 

S(t) 

position  of  liquid-solid  interface 

eff 

effective 

*,y 

Cartesian  coordinates 

1.  Introduction 

Due  to  the  attractive  features  of  latent  heat  storage,  phase 
change  materials  (PCMs)  are  mainly  used  to  store  energy  at  a  fixed 
temperature  (melting/solidification  temperature)  with  high 
energy  density  [1],  PCMs  have  been  applied  in  various  systems 
and  aspects  such  as  energy  storage  systems,  free  ventilation,  free 
heating  and  cooling  for  buildings,  spacecraft,  food,  medicine 
conservation,  smoothing  the  peaks  of  exothermic  temperature  in 
chemical  reactions  etc.  PCMs  are  a  better  choice  comparing  to  the 
sensible  heat  storage  in  applications,  because  of  its  nearly  iso¬ 
thermal  storing  mechanism  and  high  storage  density.  PCMs  absorb 
abundant  energy  through  the  phase  transition  and  release  the 
stored  energy.  Table  1  compares  the  typical  properties  of  different 
thermal  storage  materials  tested  at  the  room  temperature. 
It  indicates  that  PCMs  can  save  92.8%  of  mass  and  up  to  90%  space 
to  store  the  same  amount  of  thermal  energy  comparing  to  the 
sensible  thermal  materials  such  like  concrete  and  water. 

Therefore,  PCMs  attract  the  researchers'  interesting  in  practically 
individual  and  incorporation  applications  in  different  engineering 
fields.  The  possibility,  feasibility,  thermal  performance  and  economic 


Table  1 

Properties  comparison  of  different  thermal  storage  materials. 


Property 

Unit 

Materials 

Concrete  (sand 
and  gravel) 

Water 

Organic 

PCM 

Inorganic 

PCM 

Density 

[kg/m3] 

2240 

1000 

800 

1600 

Specific  heat 
capacity 

[kj/kg  K] 

1.1 

4.2 

2.0 

2.0 

Latent  heat 

[kj/kg  K] 

- 

- 

190 

230 

Storage  mass  for 
106  kj,  avg 

[kg] 

60,000 

16,000 

5300 

4350 

Storage  volume  for 
106  kj,  avg 

[m3] 

26.8 

16 

6.6 

2.7 

Relative  storage 

mass 

13.79 

4 

1.25 

1.0 

Relative  storage 
volume 

10 

6 

2.5 

1.0 

1.  Storage  mass  and  volume  are  calculated  for  storing  106  kj  energy  with  a 
temperature  rise  of  15  K  for  concrete  (sand  and  gravel)  and  water. 

2.  Relative  mass  and  volume  are  based  on  the  inorganic  phase  change  materials. 


analysis  of  using  PCMs  call  a  series  of  theoretical  and  experiential 
investigations,  the  mainly  two  methods  used  to  research  the  thermo¬ 
physical  parameters  of  PCMs  and  interrelated  systems. 

The  experimental  approaches  offer  a  better  indication  of  the 
actual  PCM  behavior  and  performance  in  comparison  to  theore¬ 
tical  analysis,  as  the  experimental  tests  can  present  the  PCMs’ 
behaviors  more  directly,  visibly  and  credibly  without  any  pre-set 
assumptions  however,  the  experiments  are  unachievable  in  some 
cases,  such  as  the  large  scale  or  unsteady  around  environment,  so 
there  are  still  a  few  points  need  to  be  considered: 

•  Time  and  cost  consuming  will  be  higher  than  a  theoretical 
approach,  hence,  the  budget  and  processes  needs  to  be  estab¬ 
lished  and  scheduled. 

•  The  scales  of  the  experiment  test  rig  have  to  be  considered,  and 
then  suitable  laboratory  location  and  space  needs  to  be 
determined  to  house  the  rig. 

•  Test  rig  needs  to  be  constructed  and  operated  properly. 

•  Proper  environment  parameters  have  to  be  simulated  and 
adjusted  to  imitate  the  practical  environment. 

•  Relevant  parameters  need  to  be  monitored,  measuring  appara¬ 
tuses  need  to  be  calibrated,  and  the  failed  data  need  to  be 
eliminated. 

Except  these  agonizing  experimental  matters,  there  are  still 
some  unavoidable  testing  errors.  However  the  theoretical  methods 
can  avoid  all  these  weakness  and  predicate  the  PCMs  performance 
conveniently.  The  major  advantage  of  the  theoretical/numerical 
approaches  is  that  various  conditions  can  be  carried  out  by 
changing  the  variables  in  a  numerical  model.  Therefore,  more 
and  more  investigators  prefer  to  study  the  phase  change  problems 
by  mathematical  solutions  and  numerical  simulations.  There  are 
only  eight  review  papers  on  PCMs  that  involves  the  aspect  of 
mathematical  solutions  and/or  numerical  modeling  of  latent  heat 
thermal  energy  storage  (LHTES)  have  been  published  during 
2002-2012.  Table  2  summaries  these  eight  relative  review  papers 
and  comments  their  mathematical  and/or  numerical  aspects  of 
LHTES.  However,  few  logistic  and  comprehensive  reviews  on  the 
mathematical  solutions  and  numerical  modeling  for  melting  and 
solidification  processes  of  PCMs  were  found  in  published  litera¬ 
tures  even  though  there  are  many  individual  publications  in  phase 
change  problems,  particularly  in  heat  transfer  mechanisms  and 
simulations.  Due  to  the  complexity  of  phase  transformations, 
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Table  2 

Some  review  papers  related  to  numerical  study  of  phase  change  problems. 


Reference 

Published 

year 

Journal 

Comment 

Zalba  et  al.  [1] 

2003 

Applied  thermal 
engineering 

The  review  was  divided  into  three  parts:  materials,  heat  transfer  and  applications.  Heat  transfer 
was  considered  mainly  from  a  theoretical  point  of  view,  considering  different  simulation  techniques. 

Verma  et  al.  [2] 

2008 

Renewable  and  sustainable 
energy  reviews 

Mathematical  models  of  a  LHTES  were  reviewed  to  optimize  the  material  selection  and  to  assist 
in  the  optimal  designing  of  the  systems.  Some  important  characteristics  of  different  models  and  their 
assumptions  used  were  presented  and  discussed  as  well. 

Jegadheeswaran 
and  Pohekar  [3] 

2009 

Renewable  and  sustainable 
energy  reviews 

Various  thermal  conductivity  enhancement  techniques  reviewed  in  this  paper,  and  the  issues  related 
to  mathematical  modeling  of  LHTES  with  enhancement  techniques  are  also  discussed. 

Zhu  et  al.  [4] 

2009 

Energy  conversion  and 
management 

This  paper  presented  an  overview  on  dynamic  characteristics  and  energy  performance  of  buildings 
employing  PCMs  by  three  research  methods  used,  i.e.,  simulation,  experiment,  combined 
simulation  and  experiment 

Agyenim  et  al.  [5] 

2010 

Renewable  and  sustainable 
energy  reviews 

This  paper  provided  the  formulation  of  the  phase  change  problem.  In  terms  of  problem  formulation, 
it  was  concluded  that  the  common  approach  has  been  the  use  of  enthalpy  formulation. 

Kuznik  [6] 

2011 

Renewable  and  sustainable 
energy  reviews 

The  review  was  the  first  comprehensive  review  of  the  integration  of  phase  change  materials  in  building 
walls.  Various  numerical  studies  concerning  the  integration  were  summarized. 

Dutil  et  al.  [7] 

2011 

Renewable  and  sustainable 
energy  reviews 

The  review  presented  models  based  on  the  first  law  and  on  the  second  law  of  thermodynamics.  This  paper 
tried  to  enable  one  to  start  his/her  research  with  an  exhaustive  overview  of  the  subject. 

Zhou  [8] 

2012 

Applied  energy 

This  paper  reviewed  previous  works  on  latent  thermal  energy  storage  in  building  applications,  including 
PCMs,  current  building  applications  and  their  thermal  performance  analyses,  as  well  as  numerical 
simulation  of  buildings  with  PCMs. 

a  good  understanding  of  the  heat  transfer  mechanisms,  phase 
change  characteristics  and  the  differences  among  various  mathe¬ 
matical  and  numerical  simulation  methods  is  required  before 
researchers  starting  their  theoretical  study.  The  aim  of  this  paper 
is  to  comprehensively  study  the  mathematical  and  simulation 
modeling  applied  for  LHTES.  Firstly  it  discusses  the  issues  of  the 
heat  transfer  mechanisms  during  the  charging  and  discharging 
processes  and  the  Stephan  problem  and  Neumann  problem. 
Secondly,  it  sums  the  strong  and  weak  aspects  of  various  mathe¬ 
matical  solutions  when  conduction  and  convection  heat  transfer 
occurring  inside  the  PCMs  by  the  considering  fixed  grid  method 
and  the  adaptive  grid  method.  Finally,  numerical  modeling  of 
various  PCM  packages'  geometries  in  terms  of  macro-encapsulated 
PCMs  and  micro-encapsulated  PCMs  are  studied. 


2.  Heat  transfer  mechanisms  during  the  phase 
transformations 

2.1.  Conduction  and  convection  heat  transfer 

During  the  charging  and  discharging  processes,  the  possible 
heat  transfer  mechanisms  are  conduction  and  convection.  How¬ 
ever,  the  issue  of  which  heat  transfers  mechanism  takes  the  main 
role  in  different  stages  of  phase  transformation  has  been  argued 
for  decades.  When  PCMs  are  used  to  store  or  release  thermal 
energy,  conduction  is  usually  believed  playing  the  most  important 
role  on  the  heat  transfer  during  the  phase  transformation  pro¬ 
cesses  [2-4].  As  early  as  1831,  Lame  and  Clapeyron  have  conducted 
the  first  study  on  phase  change  problem  by  only  considering  the 
pure  conduction.  However,  some  researchers  persist  that  natural 
convection  is  a  more  important  mechanism  in  the  phase  change 
process  especially  in  the  melt  region.  Lazaridis  [100]  proposed  a 
study  of  the  relative  importance  of  conduction  and  convection  in 
1970.  A  pioneer  study  performed  by  Sparrow  et  al.  in  1977  [41], 
they  concluded  in  their  study  that  natural  convection  could  not 
be  ignored  in  the  analysis  of  phase  change  problems.  In  1994, 
Hasan  [5]  concluded  that  the  convection  heat  transfer  active  an 
important  role  in  the  melting  process,  and  a  simplified  model  by 
only  considering  the  conduction  heat  transfer  does  not  describe 
the  melting  process  properly.  Later,  Lacroix  et  al.  [6]  reported 
the  similar  findings  in  their  research  that  natural  convection  is 
the  main  heat  transfer  mechanism  during  the  melting  process. 


In  1999,  Velraj  [7]  obtained  the  same  conclusion  in  their  research 
and  reported  that  during  the  melting  process  natural  convection 
occurs  in  the  melt  layer  which  results  in  the  heat  transfer  rate 
increase  comparing  to  the  solidification  process.  Buddhi  et  al.  [8] 
proposed  an  explanation  for  this  phenomenon  that  the  density 
differences  between  the  solid  and  liquid  PCM  induce  the  buoy¬ 
ancy,  which  causes  the  convective  motions  in  the  liquid  phase. 
However,  Zhang  and  Yi  [9]  believed  that  with  the  solid  PCM 
melting  into  liquid,  the  PCM  volume  keeps  increasing,  which 
results  in  the  convection  becoming  the  predominant  heat  transfer 
mechanism  rather  than  conduction.  Sari  et  al.  [10]  found  that  the 
heat  transferred  from  a  heat  exchanger  to  a  PCM  of  stearic  acid  is 
largely  influenced  by  natural  convection  at  the  melting  layer,  in 
addition  to  conduction  and  forced  convection  heat  transfer. 

For  the  melting  process,  the  PCM  changes  its  phase  from  solid 
to  a  mushy  state,  and  then  liquid,  which  is  reversible  during  the 
solidification  process.  Hence  there  are  six  stages  to  finish  a 
charging  and  discharging  cycle.  Therefore  it  is  possible  in  certain 
stages  of  the  phase  transformation  process  there  are  more  than 
one  kind  of  heat  transfer  mechanisms  acting  at  work,  but  how  to 
weight  the  percentages  of  conduction  and  convection  heat  transfer 
in  each  stage  has  been  the  main  challenge  for  the  researchers 
currently.  This  paper  will  review  the  most  common  research 
methods  used  for  the  conduction  and  convection  heat  transfer 
inside  the  PCM  packages. 

2.2.  Stefan  problem 

Moving  boundary  problem  named  Stefan  problem  is  another 
issue  to  develop  a  numerical  modeling  of  PCMs  [11,12],  The 
simplest  of  the  Stefan  problems  is  the  one-phase  Stefan  problem 
since  only  one-phase  involved.  The  term  of  ‘one-phase’  designates 
only  the  liquid  phases  active  in  the  transformation  and  the  solid 
phase  stay  at  its  melting  temperature.  Stefan’s  solution  with 
constant  thermophysical  properties  shows  that  the  rate  of  melting 
or  solidification  in  a  semi-infinite  region  is  governed  by  a  dimen¬ 
sionless  number,  known  as  the  Stefan  number  (St), 

StJW'-TM  (1) 

where  Q  is  the  heat  capacity  of  the  liquid  PCMs,  L  is  the  latent 
heat  of  fusion,  and  T(  and  Tm  are  the  surrounding  and  melting 
temperatures,  respectively. 
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The  most  published  solutions  are  to  solve  the  one-dimensional 
cases  subjected  to  an  infinite  or  semi-infinite  region  with  simple 
initial  and  boundary  conditions.  Whereas  with  the  heat  transfer 
continuing,  the  interface  boundary  is  constantly  moving  as  the 
liquid  and  solid  phases  shrinking  and  growing,  which  disable  the 
prediction  of  the  boundary  location  [13],  Because  of  that  the  solid- 
liquid  interface  is  not  fixed,  but  moving  with  time,  the  heat 
transfer  mechanisms  during  a  PCM  phase  transformation  process 
are  complex.  Therefore  the  phase  change  transition  is  difficult  to 
analyze  owing  to  the  three  reasons:  the  solid-liquid  interface  is 
moving;  the  interface  location  is  non  linear;  it  consists  of  thermal 
conduction  and  natural  convection  heat  transfer  mechanisms. 
Since  these  three  factors,  the  non-linearity  of  the  governing 
equations  is  introduced  to  the  moving  boundary,  and  the  precise 
analytical  solutions  are  only  possible  for  a  limited  number  of 
scenarios  [14].  This  view  is  shared  by  Kurklu  et  al.  [15]  who 
pointed  out  that  the  mathematical  models  have  been  proposed  in 
each  research  only  applied  to  very  specific  boundary  conditions, 
hence  are  not  feasible  to  more  complex  practical  applications.  This 
Stefan  problem  is  further  complicated  by  the  fact  that  most  of  the 
methods  previously  proposed  by  researchers  only  involve  one 
moving  boundary,  whereas  it  actually  consist  of  more  than  one 
interface  location  [16], 


2.3.  Neumann  problem 


The  Stefan  problem  was  extended  to  the  two-phase  problem, 
the  so-called  Neumann  problem  which  is  more  realistic  [17], 
In  Neumann  problem,  the  initial  state  of  the  PCM  is  assumed  to 
be  solid,  during  the  melting  process,  its  initial  temperature  does 
not  equal  to  the  phase  change  temperature,  and  the  melting 
temperature  does  not  maintain  at  a  constant  value.  If  the  melting 
happens  in  a  semi-infinite  slab  (0<x<oo),  the  solid  PCM  is 
initially  at  a  uniform  temperature  (Ts  <  Tm),  and  a  constant 
temperature  (T0)  is  imposed  on  the  slab  surface  x  =  0,  with  the 
assumptions  of  constant  thermo-physical  properties  of  the  PCM, 
the  problem  can  be  mathematically  expressed  as  follows: 

Heat  conduction  in  solid  or  liquid  region 


pc 


dT  _  d2T 
at-  dx2 


(2) 


where  p  is  the  density,  C  is  the  specific  heat,  k  is  the  thermal 
conductivity,  and  t  and  x  are  the  time  and  space  coordinates 
respectively. 

The  heat  fluxes  transferring  from  the  liquid  phase  to  the  solid- 
liquid  interface,  as  well  as  the  latent  heat  absorbing  rate  by  the 
melting  PCM,  the  movement  of  the  solid/liquid  interface  can  be 
determined  through  the  following  energy  balance: 


dxt  dT/  dTs 
pL-dt  =  ~k^+k^ 


(3) 


where  l  is  the  latent  heat  of  fusion  of  the  PCM  and  p  is  the  density 
of  liquid  PCM. 

Initial  condition 

T(x>0,t  =  0)  =  rs<Tm  (4) 

Solid-liquid  interface  temperature 
T(x  =  S(t),  t>0)=Tm  (5) 

With  the  following  boundary  conditions: 

T( 0,  t)  =  T0>Tm  for  t>0  (6) 


T(x,  t)  =  Ts  for  x->oo,  f>0  (7) 

where  S(t)  is  the  position  of  the  solid-liquid  interface  (melting 
front).  Fig.  1  clearly  illustrates  this  two-phase  Stefan  problem. 


T=Tm 


Fig.  1.  Schematic  illustration  of  the  two-phase  Stefan  problem. 


Analytical  solution  to  such  a  problem  was  obtained  by  Neu¬ 
mann  in  term  of  a  similarity  variable  p 


(8) 


The  final  Neumann's  solution  can  be  written  as 
Interface  position 


<5(f)  =  2 X\J  oc  it 


Temperature  of  the  liquid  phase 


T(x,t)  =  T,-(T,-Tm) 


er/(x/2V  cc  it) 
erf(A) 


(9) 

(10) 


Temperature  of  the  solid  phase 


T(x,t)  =  Ts+{Tm-Ts) 


erfc(x/2s/ozst) 
erfc(As/  OC  //  OCs) 


(11) 


The  A  in  Eq.  (9)— ( 11 )  is  the  solution  to  the  transcendental 
equation 

- S.tl - -  =/U/^ 

exp  (A  )erf (A)  exp(  cc  ,A'  /  cc  s)erfc(Ay/  oc  (/  oc  s) 

(12) 

where 

Stl  =  Ci<Tt-Tm)  (13) 

Sts  =  Cs(Tm~Ts)  (14) 

However,  the  Neumann's  solution  is  applicable  only  for  moving 
boundary  problems  in  the  rectangular  coordinate  system. 

2.4.  Other  possibility  mechanisms  in  phase  change  process 

Furthermore,  there  are  several  other  issues  with  the  use  of  a 
theoretical  approach  in  the  study  of  PCMs.  Alexiades  [18]  pointed 
out  that  there  were  many  mechanisms  involved  during  a  PCM 
phase  transition,  such  like  a  change  in  volume,  density,  thermal 
conductivity,  specific  heat  capacity,  super-cooling,  etc.  Conse¬ 
quently,  accurate  reflection  of  the  proposed  mathematical  solution 
and  numerical  model  need  to  consider  the  dynamic  properties  of 
the  phase  change  process.  Another  major  issue  with  PCMs  is  that 
they  act  as  self-insulating  materials.  When  PCM  solidification 
occurs  from  the  top  of  the  heat  surface,  solid  insulating  layer  will 
be  developed  which  moves  inward  during  the  whole  solidification 
process.  With  the  increase  in  the  size  and  thickness  of  the  solid 
layer,  the  heat  transfer  rate  from  the  liquid  PCM  to  the  heat 
exchanger  surface  decreases  until  it  becomes  so  small  that  will  not 


S.  Liu  et  al.  /  Renewable  and  Sustainable  Energy  Reviews  33  (2014)  659-674 


663 


be  possible  to  maintain  at  an  acceptable  heat  transfer  rate.  The 
self-insulating  characteristic  of  the  PCM  becomes  a  strong  problem 
when  bulk  containerization  is  used.  For  example  the  large  cylind¬ 
rical  used  as  containers  filled  with  PCM,  the  self-insulating  issue 
can  become  veiy  serious.  Large  heat  exchange  area  will  accelerate 
the  formation  of  the  crust  and  heavily  decrease  the  heat  transfer 
rate.  Radhakrishnan  et  al.  [14]  analyzed  that  the  method  to  lighten 
the  crust  issue  is  minimizing  the  depth  of  the  PCM  pack.  So  that 
even  a  large  majority  of  the  PCM  performs  solidification,  the 
thermal  resistance  is  reduced  to  allow  the  removal  of  the  latent 
heat  released  from  the  solidified  PCM.  The  self-insulating  char¬ 
acteristic  of  the  PCM  during  discharging  process  needs  to  be 
considered  in  any  mathematical  model  as  this  can  affect  thermal 
energy  discharging  time  and  rate.  Furthermore,  the  self-insulating 
has  a  direct  effect  on  the  phase  transition  boundary,  where  the 
solid  and  liquid  region  meets  and  co-exists,  which  is  linked  back  to 
the  Stephan  Problem  and  Neumann's  method  again. 

Logical  heat  transfer  mechanisms  during  the  charging  and 
discharging  periods  of  the  PCM  are  the  essential  elements  to 
develop  any  accurate  mathematical  model.  However,  as  described 
there  are  quite  a  few  issues  to  properly  describe  the  heat  transfer 
mechanisms:  conduction  and  natural  convention  during  the  phase 
change  process.  This  paper  reviews  most  of  the  researches 
employing  the  mathematical  calculation  and  numerical  simulation 
published  in  the  articles  and  reports,  but  does  not  find  any 
consensuses  on  which  one  is  the  dominant  heat  transfer  mechan¬ 
ism,  which  one  can  be  ignored  and  how  to  combine  these  two 
transfer  mechanisms  in  each  modeling.  Different  researchers  have 
developed  their  specific  methods  to  evaluate  the  thermal  energy 
transfer  during  the  phase  transition  process  for  various  situations 
[19-23], 


3.  Mathematical  solutions 
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Fig.  2.  Position  of  the  moving  boundary  in  a  fixed  grid. 


moving  boundary  will  be  located  between  two  adjacent  grid  points, 
for  instance,  between  mAx  and  (m  +  1) Ax  (at  a  location  of  e Ax, 
where  0  <  e  <  1),  as  illustrated  in  Fig.  2. 

The  mathematical  solution  of  the  one-phase  ice-melting  pre¬ 
sents  a  simple  illustration  of  this  method.  The  instantaneous 
temperature  at  point  (i,j+l)  is  calculated  from  the  temperatures 
of  the  previous  step  on  the  basis  of  the  following  formulation: 

Ty  +1  =  r,j  +  ( r  i  - 1  j  -  2  r,- j  +  r, + 1 J )  i  =  0,1,...m  —  1.  (15) 


In  terms  of  three-point  Lagrange  interpolation  instead  of 
Eq.  (14),  the  temperature  at  x=  mAx  is  computed  as, 


T 


m,n  + 1 


—  Tmn 


m-l,n 


(16) 


The  analyzes  of  the  heat  transfer  mechanisms  in  phase  change 
processes  are  complex  owing  to  the  movable  solid-liquid  bound¬ 
ary,  which  depends  on  the  latent  heat  absorbing  or  releasing  rate 
at  the  interface.  In  order  to  simplify  the  simulation  of  the  phase 
change  problems,  assumptions  have  been  made  by  various 
researchers  in  their  studies.  In  this  section,  different  numerical 
simulations  in  terms  of  conduction  heat  transfer  and  conduction 
and/or  convection  heat  transfer  are  discussed  and  analyzed. 

3.1.  Conduction  acting  as  the  only  heat  transfer  mechanism 

To  develop  a  mathematical  model,  the  dominant  heat  transfer 
mechanism  whether  conduction  or  convection,  during  the  char¬ 
ging  and  discharging  processes,  need  to  be  analyzed  and  deter¬ 
mined.  Although  there  are  quite  a  few  issues  in  the  heat  transfer 
mechanisms  during  the  phase  transformation,  some  researchers 
simplify  this  process  by  only  considering  conduction  in  pure  sub¬ 
stances.  In  this  case,  the  mathematical  solutions  can  be  classified 
as  follows:  fixed  grid  and  adaptive  grid. 

3.1.1.  Fixed  grid 

Fixed  grid  method  does  not  require  explicit  treatment  of 
conditions  on  the  phase  change  boundary,  and  just  requires  the 
use  of  fixed-grid  schemes  able  to  cope  with  strong  nonlinearities. 
Therefore  this  method  is  able  to  utilize  standard  solution  proce¬ 
dures  for  the  energy  equations.  They  are  amenable  to  physical 
interpretation  and  easy  to  implement,  especially  for  multidimen¬ 
sional  problems  and  for  mufti-interface  problems. 

In  the  fixed  grid  method,  the  heat  flow  equation  is  approximated 
by  finite  difference  replacements  for  the  derivatives  in  order  to 
calculate  the  values  of  temperature  T,j,  at  any  time  t„  =  nAt,  the 


The  variation  of  the  location  of  the  moving  boundary  is 


£n  + 1  —  £n  ~ 


m-  l,n 


(17) 


As  shown,  the  numerical  solution  of  this  fixed  grid  method  is 
carried  out  on  a  space  grid  remaining  at  a  fixed  value  throughout 
the  mathematical  calculation.  To  approximate  the  Stefan  condi¬ 
tions  both  for  the  moving  boundary  and  the  partial  differential 
equation  at  the  adjacent  grid  points,  various  mathematical  solu¬ 
tions  have  been  proposed.  Murray  and  Landis  [24]  proposed  two 
fictitious  temperatures:  one  is  achieved  by  quadratic  extrapolation 
from  the  temperatures  in  the  solid  PCM  and  the  other  one  is  from 
the  temperatures  in  the  liquid  PCM.  The  melting  temperature  and 
the  position  of  the  moving  interface  are  incorporated  in  the 
fictitious  temperatures,  which  are  then  substituted  into  a  stan¬ 
dard  approximate  to  calculate  the  temperature  near  the  inter¬ 
face  instead  of  using  special  formulae.  Ciment  and  Guenther  [25] 
introduced  a  method  of  spatial  mesh  refinement  on  both  sides  of 
the  moving  boundary.  With  this  method,  Lazaridis  [26]  applied  the 
explicit  finite  difference  approximations  to  solve  two-phase  soli¬ 
dification  problems  in  both  two  and  three  space  dimensions. 

The  major  advantage  of  the  fixed  grid  method  is  that  the 
mathematical  calculation  of  the  phase  transition  can  be  achieved 
through  simple  modifications  of  existing  heat  transfer  numerical 
methods  and/or  software.  As  such,  they  can  be  used  for  modeling  of 
a  variety  of  complex  moving  boundary  problems.  Basu  and  Date  [27] 
and  Voller  et  al.  [28]  reviewed  the  applications  of  the  fixed  grid 
methods  for  phase  change  problem.  The  two  mainly  fixed-grid 
methods  that  have  been  proposed  are  enthalpy-based  methods 
and  the  effective  heat  capacity  method.  However,  the  fixed  grid 
method  sometimes  cannot  work  efficiently  when  the  boundary 
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Fig.  3.  Relationships  between  enthalpy  and  temperature  for  isothermal  phase 
change. 


moves  along  a  distance  larger  than  a  space  increment  in  a  time  step, 
depending  on  the  velocity  of  the  moving  boundary.  Therefore,  the 
adaptive  grid  is  developed  by  researchers  and  will  be  discussed  in 

Section  3.1.2. 


3.I.I.I.  Enthalpy  based  method.  The  existence  of  nonlinearity  of  the 
liquid-solid  interface  and  the  unknown  location  of  the  moving 
boundary  are  the  two  main  challenges  in  simulating  the  phase 
change  problems  [291.  To  deal  with  these,  a  new  model  formu¬ 
lation  called  the  enthalpy  method  was  introduced  [30].  In  the 
enthalpy  models,  the  enthalpy  is  used  as  a  dependent  variable 
value  along  with  temperature.  By  introducing  an  enthalpy  method, 
the  phase-change  problem  becomes  much  easier  since  the 
governing  equations  are  the  same  for  the  two  phases.  The 
interface  conditions  are  automatically  achieved  and  a  mushy 
zone  between  the  two  phases  is  created.  This  zone  avoids  sharp 
discontinuities  that  may  create  some  numerical  instability.  Hence, 
the  enthalpy  mathematical  methods  are  most  attractive  and 
commonly  used  to  handle  phase  change  problems  where  a 
solid-liquid  mush  zone  is  present  between  two  phases.  Hunter 
[31]  confirmed  that  the  enthalpy  method  is  most  suitable  for 
typical  applications,  the  reason  for  that  is  this  method  does  not 
require  explicit  treatment  of  the  conditions  on  the  phase  change 
boundary  [32], 

Enthalpy-based  method  can  be  illustrated  by  considering  a 
one-dimensional  heat  conduction  controlled  phase  transforma¬ 
tion.  For  a  phase  change  process,  energy  conservation  can  be 
expressed  in  terms  of  total  volumetric  enthalpy  and  temperature. 
This  relationship  between  the  enthalpy  and  temperature  is  usually 
assumed  to  be  a  step  function  for  isothermal  phase  change  problems. 
Fig.  3  shows  the  enthalpy-temperature  curves  for  isothermal  phase 
change. 

For  isothermal  phase  change,  the  conservation  of  energy  can  be 
expressed  in  terms  of  temperature  and  the  total  volumetric 
enthalpy  [33] 

^  =  V[lc(VT)]  (18) 

ot 

where  the  total  volumetric  enthalpy  H  is  the  sum  of  sensible  and 
latent  heat,  and  can  be  expressed  in  terms  of  temperature  of  the 
PCM  as  follows: 

H(T)  =  AH(T)+Pif(T)A  (19) 

The  first  term  on  the  right  side  of  Eq.  (19)  accounts  for  the  sensible 
heat  of  the  PCM  while  the  second  term  accounts  for  the  latent  heat 
of  the  PCM. 


And  where 
AH  =  f  pcdT 


(20) 


To  be  able  to  calculate  the  total  enthalpy,  the  liquid  fraction  A  is 
given  as  follows: 


!0  T  <Tm  solid 
]0, 1[  T  =  Tm  mushy 
1  T  >Tm  liquid 


(21) 


Integrating  the  Eqs.  (18)  and  (20),  an  alternate  form  for  one¬ 
dimensional  heat  transfer  in  the  PCM  was  obtained: 


dAH 

dt 


d  (adAH) 
dx\  dx  J 


(22) 


where  a  is  the  thermal  diffusivity. 


3.3. 3.2.  Effective  heat  capacity  method  ( energy  based  method).  The 
effective  heat  capacity  method  is  also  a  common  method  used  in 
modeling  phase  change  problems  except  for  the  enthalpy  method. 
The  effective  heat  capacity  method  was  first  proposed  in  [34,35],  In 
the  effective  heat  capacity  method,  the  latent  heat  is  approximated 
by  a  large  heat  insensible  form  over  the  phase  change  temperature 
interval,  (T2  — T i)  [36j.  The  effective  heat  capacity  of  the  PCM  ( ceff ) 
is  directly  proportional  to  the  energy  stored  and  released  during  the 
phase  change  but  inversely  proportional  to  the  interval  of  the 
melting  or  solidification  temperature  range.  During  the  phase 
change  the  heat  capacity  of  the  PCM  is 

Ceff  =  T  L  T  +Cs  (23) 

1  2~  1  1 

where  T]  is  the  temperature  at  which  melting  or  solidification 
begins  and  T2  is  the  temperature  at  which  the  PCM  is  totally  melted 
or  solidified.  The  following  is  a  definition  of  the  effective  heat 
capacity  for  each  phase  change  period. 

(cs  T<Ti  solid 

tj^+Cs  T,  <T<T2  mushy  (24) 

C;  T>T2  liquid 


In  terms  of  the  definition  of  the  effective  heat  capacity,  the 
energy  equation  for  one  dimension  can  be  written  as  follows: 


dT  ,  d2T 
pCe"Tt=  kW 


(25) 


Although  the  enthalpy  based  method  and  the  effective  heat 
capacity  method  achieve  much  simpler  mathematical  operations 
with  reasonable  accuracy,  it  is  worth  noting  that  each  of  the 
methods  has  its  inherent  disadvantages.  Enthalpy  method  is 
difficult  to  handle  super-cooling  and  temperature  oscillation 
problems.  Whilst  the  effective  heat  capacity  method  is  difficult 
in  resolving  the  phase  change  problems  when  the  phase  change 
temperature  range  is  small,  and  is  not  applicable  for  the  cases 
where  phase  change  occurs  at  fixed  temperature. 

Two  approaches,  finite-difference  and  finite-element  techniques, 
are  firstly  used  to  solve  the  phase  change  problems  numerically. 
The  first  finite-difference  studies  were  carried  out  by  Lame  and 
Clapeyron  in  1831  and  Stefan  in  1891,  concerning  ice  formation. 
In  1912  the  results  of  Neuman  were  published,  establishing  the  precise 
solutions  to  more  general  phase  change  problems.  In  1943,  London 
and  Seban  [37]  analyzed  the  process  of  ice  formation  for  different 
geometries  (cylinder,  sphere,  and  flat  plate).  However,  Shamsundar 
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and  Srinivasan  [38]  asserted  that  London's  one-dimensional  formula¬ 
tion  led  to  errors  by  increasing  the  progress  of  the  solidification 
process  and  proposed  a  two-dimensional  formulation  for  cylinders.  A 
one-dimensional  formulation  was  presented  by  Bonacina  et  al.  to 
assess  the  accuracy  of  the  effective  heat  capacity  method  [39],  Rabin 
and  Korin  presented  a  multidimensional  numerical  solution  based  on 
the  effective  heat  capacity  method  for  solving  transient  phase  change 
problems  by  using  the  finite-difference  method  [40],  The  proposed 
mathematical  method  showed  a  good  agreement  with  two  exact 
solutions  and  two  numerical  solutions  based  on  finite-difference  and 
finite-element  methods. 

Although  the  finite-difference  method  has  been  traditionally 
employed  for  solving  phase  change  problems  in  the  past,  oscilla¬ 
tions  in  the  front  position  and/or  temperature  occur  due  to  a  poor 
heat  balance  when  the  phase  change  front  lies  between  nodes  in  a 
finite  difference  solution  [41,42].  Nowadays  the  finite-element 
method  has  been  more  attractive,  because  of  the  ability  of  the 
finite-element  method  to  handle  complex  coupled  thermomechanical 
mechanisms  with  various  and  complex  boundary  conditions  [43], 
Consequently,  finite-element  methods  can  either  accurately  repre¬ 
sent  the  temperature  history  or  track  the  phase  front,  or  both.  A 
few  earlier  studies  based  on  finite-element  formulations  were 
carried  out  by  Rubinsky  and  Cravahlo  in  1981  [44]  and  Ettouney 
and  Brown  in  1982  [45],  concerning  on  the  location  of  the  solid- 
liquid  interface.  To  deal  with  oscillations,  Comini  et  al.  [46]  and 
McAdie  et  al.  [47]  developed  finite-element  based  enthalpy 
formulations,  and  applied  the  three  level  Predictor-Corrector 
and  Regula-Falsi  methods  to  account  for  the  sudden  jumps  or 
sharp  changes  in  the  enthalpy-temperature  relationships.  Bhatta- 
charya  et  al.  developed  a  new  enthalpy  based  method  to  study 
phase  change  in  a  multicomponent  material  using  the  Galerkin 
finite-element  method.  This  method  can  efficiently  capture  multi¬ 
ple  thawing  fronts  which  may  originate  at  any  spatial  location 
within  the  sample  [48].  Hence,  finite-element  methods  are  pre¬ 
ferred  to  be  used  in  enthalpy  approaches  and  the  effective  heat 
capacity  method. 

Beside  the  finite-difference  and  finite-element  methods,  the 
finite-volume  and  control-volume-  finite-element  methods  were 
also  employed  to  handle  the  phase  change  problems  [39], 

3.1.2.  Adaptive  grid 

The  problem  associated  with  the  fixed  grid  method  can  be  deal 
with  by  utilizing  the  adaptive  grid  method  to  improve  the 
computational  efficiency,  using  the  so-called  adaptive  grid  (or 
front  tracking)  numerical  schemes  to  capture  the  moving  bound¬ 
ary  [49].  In  the  adaptive  grid  method,  the  exact  location  of  the 
moving  boundary  is  evaluated  on  a  grid  at  each  step.  There  are  two 
main  approaches  used  to  achieve  this:  the  interface-fitting  grids 
and  the  variable-space  grids.  The  interface-fitting  grids  (also 
referred  to  as  the  variable  time  step  methods),  where  a  uniform 
spatial  grid  but  a  non-uniform  time  step  are  used,  has  been 
repeatedly  employed  to  solve  two-phase  and  one-dimensional 
problems.  The  another  approach  is  variable-space  grids,  also 
known  as  the  dynamic  grids,  where  the  number  of  spatial  intervals 
is  kept  constant  and  the  spatial  intervals  are  adjusted  in  such  a 
manner  so  that  the  moving  boundary  lies  on  a  particular  grid 
point.  Thus,  in  these  methods  the  spatial  intervals  are  a  function  of 
time.  The  applications  of  the  variable  grid  method  to  study  the 
phase  change  problem  can  be  found  in  [50-54], 

Lacroix  and  Voller  [55]  performed  a  study  on  simulation  of 
diffusion/convection  controlled  solidification  processes  in  a  rec¬ 
tangular  cavity  by  using  the  finite-difference  techniques.  They 
concluded  that  the  fixed  grid  must  be  finer  for  solving  the  melting 
problem  of  a  material  with  a  unique  melting  temperature.  Com¬ 
parison  between  the  fixed  and  adaptive  grid  has  also  been  done  by 
Bertrand  et  al.  [56],  They  found  that  the  adaptive  grid  method  is 


better  adapted  to  the  melting  problem  than  the  fixed  grid  method. 
However  the  adaptive  grid  method  may  fail  to  simulate  the 
situations  where  the  transition  from  the  liquid  to  the  solid  phase 
is  not  a  macroscopic  surface. 

Table  3  presents  some  references  of  the  mathematical  and 
numerical  studies  undertaken  with  the  only  consideration  of  the 
conduction  heat  transfer. 

3.2.  Conduction  and  convection  heat  transfer  mechanisms 

During  the  phase  transformation  especially  melting  process, 
the  temperature  and  concentration  gradients  in  the  liquid  phase  of 
PCM  keep  varying,  which  results  in  the  movement  of  the  liquid 
PCM,  named  convection  occurring  under  the  action  of  buoyancy 
forces  due  to  the  density  gradients.  In  the  previous  numerical 
investigations,  the  convection  heat  flow  in  the  liquid  phase 
received  less  attention  than  conduction  owing  to  the  limited 
computational  capabilities  of  computer  and  the  mathematical 
complexities  to  formulate  the  convection  heat  transfer  during 
the  phase  transformation.  The  influence  of  the  natural  convection 
on  the  phase  change  problems  were  initially  considered  in  1977. 
The  study  of  the  possible  role  of  the  natural  convection  in  the  melt 
region  has  been  conducted  by  Sparrow  et  al.  in  1978  [67],  The 
findings  of  this  study  indicated  that  the  natural  convection  is  first- 
order  important  and  has  to  be  accounted  in  the  phase  transforma¬ 
tion  analysis.  A  number  of  researchers  [68-71]  have  also  reported 
that  the  convection  affects  not  only  the  rate  of  melting  or 
solidification  but  also  the  resulting  of  the  distribution  of  the  liquid 
phase  of  the  PCM  in  a  multi-component  system.  One  of  the 
pioneer  numerical  studies  was  carried  out  by  Yao  and  Chen  in 
1980  [72],  by  using  an  approximate  solution.  They  studied  the 
effect  of  the  natural  convection  on  the  phase  change  process  and 
concluded  that  it  strongly  depends  on  the  Rayleigh  number.  The 
equation  of  ke/lq  =  cRan  includes  an  effective  of  thermal  conduc¬ 
tivity  was  used  to  describe  the  influence  of  natural  convection  on 
the  melting  process  [73,74], 

The  influence  of  the  convection  heat  transfer  on  the  melting 
and  freezing  processes  has  been  studied  by  Lamberg  et  al.  [64], 
they  found  that  when  the  effect  of  natural  convection  is  omitted  in 
the  simulation.  The  numerical  result  was  poor  compared  with  the 
experimental  result  during  the  melting  process,  whilst  it  showed  a 
good  estimation  during  the  freezing  process.  Hence  it  is  essential 
to  model  natural  convection  in  the  liquid  PCM  during  the  melting 
process.  Besides  the  numerical  methods  mentioned  in  Section  3.1, 
the  integral  method  [75],  boundary  fixing  method  [76],  unstruc¬ 
tured  finite-element  method  [77],  enthalpy-porosity  approach 
[78-80],  coordinate  transformation  method  [81,82]  and  equivalent 
thermal  capacity  method  [83]  have  been  developed  for  natural 
convection  simulation. 

Yao  and  Cherney  [84]  studied  the  effect  of  the  natural  convec¬ 
tion  on  the  melting  of  a  solid  PCM  around  a  hot  horizontal  cylinder 
by  using  the  integral  method.  The  results  demonstrated  that  the 
integral  solution  had  surprising  accuracy  when  it  was  compared 
with  the  quasi-steady  solution  when  Stefan  number  was  small. 
Rieger  et  al.  [85]  investigated  numerically  the  melting  process  of  a 
PCM  inside  a  heated  horizontal  circular  cylinder.  Both  heat  con¬ 
duction  and  convection  have  been  taken  into  account  to  treat  this 
moving  boundary  problem,  and  the  complex  structure  of  the  time- 
wise  changing  physical  domain  (melt  region)  have  been  success¬ 
fully  overcome  by  applying  body  fitted  coordinate  technique. 
Ismail  and  Silva  [86]  developed  a  2D  mathematical  model  to  study 
the  melting  of  a  PCM  around  a  horizontal  circular  cylinder 
considering  the  presence  of  the  natural  convection  in  the  melt 
region.  A  coordinate  transformation  technique  was  used  to  fix  the 
moving  front.  The  numerical  predictions  were  compared  with 
available  results  to  establish  the  validity  of  the  model,  and  a 
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Table  3 

Numerical  studies  by  only  considering  the  conduction  heat  transfer. 


Reference(s)  Dimensional  Numerical  Assumptions 

simulation  method 


Comment 


Validation 


Silva  et  al.  1 D 

[57] 


Kiirklii  et  al.  2D 
[15] 


Costa  et  al.  2D 

[3] 


Gong  et  al.  2D 
[59] 


Sharma  et  al.  2D 

[60] 


Dwarka  and  3D 
Kim  [61] 


Dolado  et  al.  1 D 

[62] 


Fukai  et  al.  2D 

[63] 


Lamberg  et  al.  2D 

[64] 


Zivkovic  and  ID 
Fujii  [65] 


Enthalpy  based 
method 


Energy  based 
equation 


Energy  based 
equation. 

Fully  implicit 

finite-difference 

method. 

Enthalpy  based 
method. 

Standard  Galerkin 
finite-element 


Enthalpy  based 
method 
Fully  implicit 
finite-difference 
method 


Implicit  enthalpy 
based  method. 


Finite-difference 

method 


Enthalpy  based 
method 

Control-volume 

method 


Enthalpy  based 
method  and  Effective 
heat  capacity 
method 


Enthalpy  based 
method 


•  1-D  conduction  controlled  the  phase  The  model  can  be  used  to  predict  the  dynamic  Yes 

change  process.  performance  of  this  kind  of  vertical  rectangular 

•  The  fluid  flow  was  fully  developed  with  heat  exchange  unit  with  reasonable  accuracy 
temperature  independent  properties. 

•  Thermal  radiation  between  the  wall 
and  fluid  was  neglected. 


•  The  phases  were  homogeneous.  Calculating  the  temperature  of  air  and  PCM  Yes 

•  Conduction  heat  transfer  was  at  the  formerly  defined  volume  nodes, 

only  considered. 

•  Heat  loss  or  gain  from  the  store  was 
neglected. 


•  The  thermophysical  properties  of  the  PCM  Calculations  have  been  made  for  the  melt  fraction  Numerical 

and  fin  were  independent  of  temperature.  and  energy  is  stored  by  conduction  only  [33,58] 

•  The  PCM  was  initially  in  the  solid  phase. 

•  The  PCM  was  homogenous  and  isotropic. 


•  The  heat-transfer  fluid  was  incompressible 
and  viscous  dissipation  was  negligible. 

•  The  fluid  flow  was  radially  uniform  and 
the  axial  velocity  was  an  independent 
parameter. 

•  Thermal  losses  through  the  outer  wall 
of  the  PCM  were  negligible. 

•  Conduction  heat  transfer  was  the 
mainly  transfer  method. 

•  The  thermo-physical  properties  of  PCM's  and 
fin  material  were  independent  of 
temperature  but  different  for  solid  and  liquid 
phases. 

•  The  PCM  was  initially  in  solid  phase.  The 
PCM  was  homogeneous  and  isotropic. 

•  The  mode  of  heat  transfer  was 
conduction  only. 

•  Air  temperature  was  constant. 

•  Initial  temperatures  were  the  same  through 
the  board. 

•  There  was  no  energy  loss  to  surroundings. 
Conduction  heat  transfer  was  only 
considered. 

•  All  thermo-physical  properties  were 
constant  except  for  the  heat  capacity. 

•  Only  conduction  heat  transfer  was 
considered  inside  the  PCM  plate. 

•  The  temperature  of  the  airflow  was  analyzed 
in  a  one  dimensional  way. 

•  The  surfaces  of  the  heat  exchanger  were 
adiabatic. 

•  The  cross-sections  of  the  tubes  was  square 

•  The  convective  heat  transfer  term  was 
neglected. 

•  The  heat  conductivity  and  density  of  the  PCM 
and  container  were  constant. 

•  The  physical  properties  chosen  for  the  PCM 
were  the  average  values  of  solid  and 
liquid  PCM. 

•  The  convection  heat  transfer  was  negligible 
during  solidification  process. 

•  The  effects  of  natural  convection  within  the 
melt  were  negligible  and  can  be  ignored. 

•  The  PCM  behaved  ideally. 

•  The  PCM  was  assumed  to  have  a  definite 
melting  temperature. 


For  model,  simulation  has  been  done  by  No 

introducing 

hot  and  cold  fluids  from  the  same  end  of  the  tube. 
While  for  mode  2,  simulation  has  been  done  by 
introducing  hot  and  cold  fluids  from  the  different 
ends  of  the  tube. 


Calculations  were  made  to  study  the  effect  of  No 
thermal  conductivity  and  thermal  capacity  of  fatty 
acids,  and  conductivity  of  heat  exchanger  materials 
and  their  effect  on  melt  fraction. 


Comparison  of  the  thermal  performance  of  No 

randomly  mixed  and  laminated  PCM  drywall 
system 


A  model  was  developed  to  simulate  the  Yes 

performance 

of  a  LHTES,  and  analyze  the  heat  transfer  between 
the  air  and  a  commercial  macroencapsulated  PCM. 

A  numerical  model  predicted  well  the  experimental  Yes 
outlet  fluid  temperatures  and  the  local 
temperatures 
in  the  composite. 


Both  numerical  methods  gave  good  estimations  for  Yes 
the  temperature  distribution.  Particularly,  when 
the  temperature  range  was  2  °C,  the  effective  heat 
capacity  method  was  the  most  precise  numerical 
method. 


The  cylindrical  container  required  nearly  twice  of  Yes 
the  melting  time  as  for  the  rectangular  container  of 
the  same  volume  and  heat  transfer  area 
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Table  3  ( continued ) 


Reference(s)  Dimensional  Numerical  Assumptions 

simulation  method 


Comment 


Validation 


Saman  et  al.  2D  Enthalpy  based  •  The  thermo  physical  properties  of  PCM  The  numerical  model  was  able  to  predict  quite  Yes 

[66]  equation  can  be  assumed  constant.  accurately  the  melting  time  and  the  heat  transfer 

•  The  PCM  was  initially  solid  and  its  rate  during  the  melting  with  the  geometry  being 

temperature  was  assumed  at  a  certain  value  considered 

below  the  melting  point. 

•  The  natural  convection  was  taken 
into  account  during  charging  period. 


satisfactory  agreement  was  found.  They  concluded  that  the 
numerical  model  was  adequate  to  represent  the  physical  situation 
of  the  proposed  system. 

Furzerland  [87]  investigated  the  enthalpy  method  and  the 
coordinate  transformation  method  through  the  comparison  of 
the  solutions  of  specific  problems  of  one  dimensional  heat  transfer 
by  considering  pure  convection.  One  of  his  conclusions  was  that 
the  enthalpy  method  is  very  attractive  owing  to  these:  it  is  easy  to 
program  and  there  are  no  computational  overheads  associated 
with  tracking  the  moving  interface  within  a  specific  range  of 
fusion  temperatures. 

A  summary  for  various  numerical  solutions  on  phase  change 
problem  considering  both  conduction  and  convection  is  presented 

in  Table  4. 

4.  Numerical  models  for  various  PCMs  capsules  and  packages 

It  can  be  easily  found  that  solid-liquid  PCMs  have  been  widely 
used  for  studies  on  phase  change  problems  compared  to  solid- 
solid,  liquid-gas  and  solid-gas  PCMs.  This  is  due  to  their  large  heat 
storage  capacity  and  little  volume  change  during  the  phase  change 
process.  Since  solid-liquid  PCMs  experience  a  phase  transition 
during  the  charging  process  and/or  charging  process,  encapsula¬ 
tion  is  required  for  holding  the  liquid  and/or  solid  phases  of  the 
PCM  and  keeping  it  isolate  from  the  surrounding.  This  ensures  the 
flexibility  in  frequent  phase  change  processes,  and  an  increase  in 
heat  transfer  rate.  It  is  worth  noting  that  the  heat  transfer 
characteristics  in  the  encapsulated  PCM  would  change  signifi¬ 
cantly  depending  on  the  parameters  of  various  encapsulations. 
Hence,  vast  of  numerical  studies  on  phase  change  problems  have 
extended  to  engineering  applications.  In  this  section,  a  compre¬ 
hensive  review  on  numerical  models  for  encapsulated  PCMs  in 
terms  of  macro-encapsulated  PCM  and  microencapsulated  PCM  is 
presented. 

4.1.  Macro-encapsulated  PCMs  numerical  models 

Macro-encapsulation  is  a  common  form  of  encapsulation  used 
for  thermal  storage  applications.  The  shape  of  macro-encapsulate 
varies  from  rectangular  panels,  spheres  to  cylinders. 

4.1.1.  Rectangular  encapsulation 

The  most  intensely  investigated  LHS  containment  is  the  rec¬ 
tangular  encapsulation  because  to  its  simple  boundary  conditions 
and  easy  work  principle.  The  first  numerical  study  on  rectangular 
geometry  has  been  conducted  by  Shamsundar  and  Sparrow 
[107,108]  in  1975  and  1976  by  applying  the  finite-difference 
method  to  investigate  the  solidification  of  a  flat  plate.  Shamsundar 
also  used  the  same  method  to  the  case  of  a  square  geometry  [109] 
in  1978. 

Hamdan  and  Elwerr  [110]  presented  a  simple  2-D  numerical 
model  to  study  the  melting  process  of  a  solid  PCM  within  a 


rectangular  enclosure.  In  this  model,  the  rate  of  melting  depends 
essentially  on  the  properties  of  the  PCM,  such  as  the  thermal 
diffusivity,  viscosity,  conductivity,  latent  heat  of  fusion  and  specific 
heat.  Natural  convection  was  treated  as  the  dominant  heat  transfer 
mode  within  the  melt  region.  Conduction  was  only  assumed  to 
take  place  within  the  layer  very  close  to  the  solid  boundary.  The 
numerical  results  showed  a  high  agreement  with  previously 
published  experimental  results  [111,112].  Therefore,  it  can  be 
effectively  used  to  predict  the  energy  storage  performance  of  the 
PCMs  contained  within  the  rectangular  encapsulation. 

Further  research  was  carried  out  by  Lacroix  [113]  in  2001,  who 
presented  a  mathematical  model  based  on  the  energy  conserva¬ 
tion  equation  to  study  the  contact  melting  of  a  sub-cooled  PCM 
inside  a  rectangular  cavity  including  the  natural  convection  effect. 
Numerical  results  demonstrated  that  the  melted  fraction  at  the 
bottom  of  the  cavity  is  larger,  by  an  order  of  magnitude  than  that 
of  the  conduction  dominated  melting  at  the  top.  The  melting 
process  is  essentially  governed  by  the  magnitude  of  the  Stefan 
number  and  strongly  influenced  by  the  lateral  dimensions  of  the 
cavity.  Jiji  and  Gaye  [114]  also  analytically  examined  one¬ 
dimensional  solidification  and  melting  of  a  slab  with  uniform 
volumetric  energy  generation.  They  found  that  the  low  thermal 
conductivity  of  the  PCMs  present  a  significant  challenge  in  the 
design  of  PCM-based  electronic  cooling  systems. 

In  order  to  address  the  inherent  drawback  of  the  low  thermal 
conductivity  of  PCM,  various  numerical  solutions  have  been 
applied  for  PCM  with  different  thermal  conductivity  enhancers 
(TCEs).  Lacroix  and  Benmadda  [115]  numerically  analyzed  the 
melting  of  PCMs  in  a  rectangular  enclosure  with  horizontal  fins. 
A  2-D  enthalpy  model  was  used  to  solve  the  phase  change 
problem  using  a  fixed-gird  method.  The  effects  of  the  number 
and  length  of  fins  on  the  melting  rate  were  considered  in  this 
study.  It  was  concluded  that  a  few  longer  fins  were  more  effective 
in  the  melting  rate  than  a  large  number  of  shorter  fins.  The  effect 
of  vertical  fins  on  melting  rate  was  studied  by  the  same  team  [116] 
by  using  a  fixed-grid  enthalpy  model  in  1998.  It  was  found  that  the 
onset  of  natural  convection  in  the  melted  zone  was  delayed  when 
the  distance  between  the  fins  was  decreased.  According  to  their 
results,  the  optimum  fin  spacing  decreases  as  the  Rayleigh  number 
increases.  Akhilesh  et  al.  [117]  numerically  studied  a  rectangular 
module  with  vertical  fins,  which  was  heated  from  above.  The 
numerical  model  was  solved  by  using  finite-difference  method. 
Only  conduction  heat  transfer  was  included  in  this  study.  The 
analysis  showed  that  more  fins  increase  the  rate  of  energy  stored 
in  the  melting  PCM  until  a  critical  value.  There  is  no  obvious 
performance  enhancement  by  only  increasing  the  number  of  fins 
beyond  the  value.  Similar  results  were  obtained  in  another 
research  of  Gharebagi  and  Sezai  [118]  in  2005  considering  an 
enclosure  with  vertical  fins  added  to  a  horizontal  heated  wall.  The 
results  indicated  that  the  heat  transfer  rate  to  the  melting  PCM  can 
be  improved  by  adding  fins.  In  addition,  vertically  heated  walls 
with  horizontal  fins  exhibit  better  performance  than  horizontally 
heated  walls  with  vertical  fins. 
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Table  4 

Numerical  studies  by  considering  both  the  conduction  and  convection. 


Reference(s)  Dimensional  Numerical  solution  method  Assumptions  Comment  Validation 


Costa  et  al.  1 D 

[3] 


Halawa  et  al.  2D 

[88] 


Wang  et  al.  2D 
[90] 


Jana  [91]  2D 


Hamdan  and 
Al-Hinti 


[93] 


3D 


Zhao  et  al.  ID 

[95] 


Liuetal.  [97]  2D 
3D 


Shmuelietal.  2D 

[99] 


Energy  based  equation 
Fully  implicit  finite- 
difference  method 


Enthalpy  based  method 


Energy  based  method 
Finite-volume  method 


Enthalpy-porosity  based 
method 

Finite-volume  method 


Energy  based  equation 


Enthalpy  porosity  method 


Enthalpy  porosity  method 


•  The  thermophysical  properties  of  the 
PCM  and  fin  material  were  independent 
of  temperature. 

•  The  PCM  was  initially  in  the  solid  phase. 

•  The  PCM  was  homogenous  and  isotropic. 

•  For  the  PCM  melting  process,  the  PCM 
was  solid  and  its  temperature  was 
assumed  at  a  certain  value  below  the 
melting  point. 

•  For  freezing  process,  the  PCM  was 
initially  liquid  and  its  temperature  was 
assumed  at  certain  value  above  the 
melting  point. 

•  The  PCM  was  pure  and  homogeneous. 

•  The  melting  front  was  an  isothermal 
interface  in  reality. 

•  The  thermophysical  properties  of  the 
PCM  were  constant,  with  the  exception 
of  the  linear  density-temperature 
relation. 

•  The  natural  convection  flow  of  the  liquid 
phase  was  incompressible  and  laminar 
with  no  viscous  dissipation. 

•  The  density  of  PCM  was  constant  in 
this  study. 

•  The  density  was  considered  to  be 
constant  in  the  unsteady  and  convective 
terms  and  was  allowed  to  vary  only  in 
the  body-force  term, 

•  The  phase  change  occurred  at  a  single 
temperature, 

•  The  temperatures  of  liquid  and  solid  at 
the  interface  were  equal. 

•  The  solid-liquid  interface  was 
smooth,  plane. 

•  The  solid  phase  was  initially  at  the 
melting  temperature. 

•  The  effect  of  inertia  inside  the  thermal 
boundary  layers  was  negligible  relative 
to  that  of  buoyancy  and  viscosity. 

•  The  heat  loss  from  the  walls  of  the 
enclosure  was  negligible. 

•  The  PCM  material  was  pure,  homogenous 
and  isotropic. 

•  Boussinesq  approximation  was  used  in 
this  analysis. 

•  The  heat  transfer  process  inside  the 
molten  layer  was  one-dimensional,  and 
the  tangential  temperature  gradient  was 
neglected. 

•  The  thermo-physical  properties  of  PCM 
were  constant. 

•  The  process  was  quasi-steady,  and  the 
acting  forces  on  PCM  were  balanced  at 
every  moment  of  the  melting. 

•  The  two-phase  mixed  region  can  be 
described  with  the  porosity. 

•  The  effective  thermal  conductivity  of  the 
mixture  can  be  calculated  as  the  volume 
average  of  the  conductivities  of  porous 
matrix  material  and  PCM 

•  The  porous  media  and  fluid  were  not 
assumed  to  be  in  thermal  equilibrium. 

•  Both  solid  and  liquid  phases  were 
homogeneous  and  isotropic. 

•  The  melting  process  was  axisymmetric. 


Calculations  have  been  made  for  the  melt  Numerical 
fraction  and  energy  stored  under  conduction  [30,58] 
plus  convection  condition. 


Typical  characteristics  of  the  melting/freezing  Numerical  [89] 
of  PCM  slabs  were  discussed.  Considerations 
in  the  design  of  the  LHTES  were  also  given. 


Calculations  were  performed  to  study  the  No 
effect  of  thermal  conductivity  and  thermal 
capacity  of  fatty  acids  and  conductivity  of 
heat  exchanger  materials  and  their  effect 
on  melt  fraction. 


Prediction  of  solidification  and  melting 
processes  of  pure  materials 
using  moving  grids. 


Numerical  and 
Experimental 

[92] 


The  melting  process  of  a  solid  PCM  contained  Experimental 
in  a  rectangular  enclosure  heated  from  [94] 

a  vertical  side  at  a  constant  heat  flux 
was  investigated  analytically. 


The  expression  of  melting  velocity  was  Numerical  and 

obtained  by  using  the  lubrication  theory.  Experimental 

[96] 


The  effects  of  the  structural  parameters  of  the  Experimental 
porous  media  and  inlet  conditions  of  HTF  on  [98] 
the  thermal  performances  of  LHTES  were 
numerically  analyzed. 


The  findings  of  the  present  study  made  it  Experimental 
possible  to  define  the  heat  transfer  [100] 

mechanisms  from  the  tube  wall  to  the  liquid 
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Table  4  ( continued ) 


Reference(s) 

Dimensional 

Numerical  solution  method 

Assumptions 

Comment 

Validation 

•  The  molten  PCM  and  the  air  were 
incompressible  Newtonian  fluids, 
and  laminar  flow  was  assumed  in  both. 

PCM  and  then  to  the  solid  phase  at  various 
locations  and  instants. 

Ye  et  al. 

[101] 

2D 

Enthalpy  porosity 
approach 

•  Density  and  dynamic  viscosity  of  the 

PCM  depended  on  temperature. 

•  Variation  of  property  values  in  liquid  and 
solid  phase  PCM  was  imposed  using 
piecewise  linear  interpolation. 

•  Constant  thermophysical  properties 
were  specified  for  aluminum. 

The  numerical  study  provided  a  detailed 
knowledge  regarding  interface  heat  transfer 
rate  provides  a  deeper  understanding  the 
heat 

transfer  mechanisms. 

Experimental 

[102]  and 
Numerical 

[103] 

Tao  and  He 

[104] 

2D 

Enthalpy  method 

•  The  axial  heat  conduction  and  viscous 
dissipation  in  the  HTF  was  negligible. 

•  The  flow  of  HTF  was  treated  as  one 
dimensional  fluid  flow. 

•  The  PCM  region  was  treated  as  an 
axisymmetric  model. 

•  The  effect  of  natural  convection  of  PCM 
during  melting  was  taken  into  account 
with  an  effective  thermal  conductivity  of 
the  liquid  phase  of  PCM. 

The  effects  of  the  non-steady  inlet  conditions 
of  HTF  on  melting  time,  melting  fraction,  TES 
capacity,  solid-liquid  interface,  heat  flux  on 
tube  surface  and  HTF  outlet  temperature 
were  analyzed. 

Experimental 

[105] 

Adine  and 
Qarnia 
[106] 

2D 

Conservation  energy 
equations 

Control  volume  approach 

•  The  flow  of  HTF  was  dynamically 
developed. 

•  The  axial  conduction  and  viscous 
dissipation  in  the  fluid  were  negligible. 

•  The  thermophysical  properties  of  the 

HTF  and  PCMs  were  independent  of  the 
temperature. 

•  The  thermophysical  properties  of  solid 
and  liquid  phases  of  PCM  were  equal, 
except  their  thermal  conductivities. 

This  parametric  study  could  provide 
guidelines 

for  system  thermal  performance  and  design 
optimization. 

Experimental 

[105] 

Apart  from  applying  fins  to  PCMs,  the  effect  of  metal  foams  on 
heat  transfer  enhancement  in  rectangular  encapsulated  PCMs 
was  investigated  by  Tian  and  Zhao  [98]  in  2011.  The  numerical 
investigation  was  based  on  the  two-equation  non-equilibrium 
heat  transfer  model,  in  which  the  coupled  heat  conduction  and 
natural  convection  were  considered  at  the  phase  transition  and 
liquid  zones.  The  numerical  results  were  validated  by  experimen¬ 
tal  data.  The  main  findings  of  the  investigation  were  that  the  heat 
conduction  rate  increases  significantly  by  using  metal  foams,  due 
to  their  high  thermal  conductivities,  and  that  natural  convection  is 
suppressed  owing  to  the  large  flow  resistance  in  metal  foams. 
In  spite  of  this  suppression  caused  by  the  metal  foams,  the  overall 
heat  transfer  performance  is  improved  when  metal  foams  are 
embedded  into  PCM.  This  implies  that  the  enhancement  of  heat 
conduction  offsets  or  exceeds  the  natural  convection  loss.  The 
results  also  indicated  that  for  different  metal  foam  samples,  heat 
transfer  rate  can  be  further  increased  by  using  metal  foams  with 
smaller  porosities  and  bigger  pore  densities. 


4.1.2.  Cylindrical  containment  models 

Three  modes  of  cylindrical  PCM  container  configurations  are 
classified.  The  first  model  is  the  PCM  filling  the  shell  and  the  heat 
transfer  fluid  (HTF)  flowing  through  a  single  tube  (pipe  model). 
In  the  second  mode  the  PCM  fills  the  tube  and  the  HTF  flows  parallel 
to  the  tube  (tube-tube  model).  The  third  cylinder  model  is  the  PCM 
encapsulated  in  a  shell  and  HTF  flows  through  the  tube  outside  the 
shell  to  improve  the  heat  transfer  of  PCMs  (tube-shell  model). 

Ismail  and  Alves  [119]  numerically  analyzed  a  pipe  model  with 
the  PCM  encapsulated  in  a  shell  while  water  as  the  HTF  flowing 
through  the  tube.  Radial  conduction  was  assumed  to  be  dominated 
during  the  process  of  solidification  and  the  energy  equation  for  the 


PCM  was  solved  in  conjunction  with  the  equation  governing 
the  fluid  bulk  temperature  as  a  function  of  time  by  employing 
the  control-volume  method.  The  effect  of  Biot  number,  ratio  of 
diameters  and  inlet  fluid  temperature  on  the  thermal  performance 
of  the  storage  unit  were  presented.  Lacroix  [105]  developed  an 
enthalpy  based  model  to  predict  the  transient  thermal  perfor¬ 
mance  characteristics  of  a  pipe  LHTES  with  the  PCM  on  the  shell 
side  and  the  HTF  circulating  inside  the  tube.  Including  the  axial 
conduction  effect  in  the  PCM  and  employing  an  enthalpy  based 
method,  the  coupled  conservation  equations  governing  the  heat 
transfer  in  the  PCM  and  HTF  were  solved  iteratively  using  a  finite- 
volume  based  finite-difference  technique.  Numerical  results  were 
presented  for  various  thermal  and  geometric  parameters  such  as 
shell  radius,  mass  flow  rate  and  inlet  temperature  of  the  HTF. 
It  was  concluded  that  for  a  given  type  of  PCM,  these  parameters 
must  be  selected  carefully  in  order  to  optimize  the  performance  of 
the  storage  unit.  Then  the  tube-tube  model  was  theoretically 
studied  by  Esen  and  Durmus  [120]  in  1998.  A  theoretical  model 
based  on  an  enthalpy  method  was  developed  to  predict  the  effects 
of  various  thermal  and  geometric  parameters,  configurations  of 
the  cylindrical  container  on  the  whole  melting  time  for  different 
PCMs.  The  theoretical  results  showed  that  the  whole  melting  time 
of  PCM  depends  on  not  only  thermal  and  geometric  parameters  of 
the  container,  but  also  on  the  thermophysical  properties  of  the 
PCM.  It  was  found  that  a  thicker  PCM  pack  would  result  in  longer 
melting  time. 

During  2006,  Regin  et  al.  [121]  experimentally  and  numerically 
studied  the  third  model  that  the  melting  behavior  of  paraffin  wax 
encapsulated  in  a  cylindrical  capsule.  The  numerical  analysis 
carried  out  by  using  enthalpy  method  and  the  results  were  verified 
with  the  experimental  data.  The  experiments  have  been  done  by 
the  visualization  technique  without  disturbing  the  actual  process 
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of  melting.  Numerical  results  indicated  that  the  melting  process  is 
chiefly  governed  by  the  magnitude  of  the  Stefan  number,  phase 
change  temperature  range  and  the  capsule  radius.  The  analysis 
showed  that  the  agreement  between  the  analytical  and  experimental 
results  is  significantly  improved  when  the  results  are  obtained  by 
considering  the  phase  change  temperature  range  and  the  natural 
convection  in  the  liquid  phase  instead  of  only  considering  the 
conduction  heat  transfer. 


4.1.3.  Spherical  containment  models 

Generally,  the  PCM  in  the  sphere  is  subjected  to  constrained 
(the  solid  and  liquid  phases  have  the  same  density)  and  uncon¬ 
strained  melting  (the  solid  portion  drops  out  of  the  melting  layer 
due  to  its  higher  density). 

The  constrained  melting  of  the  PCM  within  a  spherical  container 
was  investigated  by  Khodadadi  and  Zhang  in  2001  [4],  The 
computations  were  based  on  an  iterative,  finite-volume  numerical 
procedure  by  using  the  primitive-dependent  variables,  whereby  the 
time-dependent  continuity,  momentum  and  energy  equations  in 
the  spherical  coordinate  system  were  solved.  With  the  buoyancy- 
driven  convection  gradually  dominated  the  heat  transfer  mechan¬ 
ism  due  to  the  growth  of  the  melt  zone,  the  solid  PCM  melted  faster 
in  comparison  to  the  bottom  zone.  When  the  buoyancy  effects 
became  remarkable,  as  many  as  three  time-dependent  re-circulat- 
ing  vortices  were  observed  from  the  numerical  simulation.  The 
computational  findings  were  verified  through  qualitative  con¬ 
strained  melting  experiments  using  a  high-Prandtl  number  wax. 

An  experimental  and  computational  investigation  directed  at 
understanding  the  role  of  buoyancy-driven  convection  during  con¬ 
strained  melting  of  the  PCM  inside  a  spherical  capsule  was  reported  by 
Tan  et  al.  in  2009  again  [122],  The  computations  were  based  on  an 
iterative,  finite-volume  numerical  procedure  that  incorporated  a 
single-domain  enthalpy  formulation  for  simulation  of  the  phase 
change  phenomenon.  A  Darcy's  Law-type  porous  media  treatment 
linked  the  effect  of  the  phase  change  on  convection.  The  computa¬ 
tional  predictions  pointed  to  the  strong  thermal  stratification  in  the 
upper  half  of  the  sphere  that  resulted  from  rising  of  the  molten  PCM 
along  the  inner  surface  of  the  sphere  thus  displaced  the  colder  fluid. 
The  measured  temperature  data  and  computational  results  near  the 
bottom  indicate  the  establishment  of  an  unstable  fluid  layer  that 
promotes  chaotic  fluctuations  and  is  responsible  for  waviness  of  the 
bottom  of  the  PCM.  Regin  et  al.  [123]  also  discussed  an  experimental 
and  numerical  study  of  a  spherical  capsule  filled  with  a  paraffin  wax 
(melting  temperature  range  of  52.9-61.6  °C)  as  the  constrained  PCM 
and  placed  in  a  convective  environment  during  the  melting  process. 
Time-resolved  temperature  readings  were  obtained  by  using  a  rake  of 
nine  thermocouples  that  were  placed  on  a  horizontal  plane  passing 
through  the  center  of  the  sphere.  Keeping  the  initial  temperature  of 
the  PCM  constant,  the  temperature  of  the  convective  heat  transfer 
fluid  varied  from  70  to  82  °C  corresponding  to  Stefan  numbers  of 
0.1143-0.2501.  A  ID  model  for  phase  change  by  employing  the 
enthalpy  method  and  finite-difference  formulation  was  developed. 
The  model  was  then  used  to  investigate  the  effect  of  capsule  size  and 
the  Stefan  number  on  the  instantaneous  heat  flux,  liquid  fraction  and 
melting  time. 

Considering  the  spherical  geometries,  apart  from  several  ana¬ 
lytical  investigations  of  diffusion  controlled  phase  change  pro¬ 
cesses  reported  before  1980s,  the  first  study  on  the  unconstrained 
melting  of  the  PCM  within  spherical  container  was  carried  out  by 
Moore  and  Bayazitoglu  [124]  in  1982.  A  mathematical  model  was 
developed  by  assuming  that  the  liquid  remains  at  the  fusion 
temperature  during  the  whole  period  and  solved  by  employing 
the  perturbation  method.  The  solid-liquid  interface  positions  and 
the  temperature  profiles  were  determined  for  various  Stefan  and 
Fourier  numbers.  The  experimental  data  agreed  well  with  the 


simulation  results.  During  1987,  Roy  and  Sengupta  [125]  per¬ 
formed  a  similar  study  and  found  an  analytical  solution  for  the 
melting  rate  at  the  lower  surface  of  the  solid  PCM  in  contact  with 
the  heated  surface.  Later,  Roy  and  Sengupta  [126]  performed 
another  study  on  the  unconstrained  melting  process  inside  a 
sphere.  Two  zones  of  the  melting  were  identified:  a  thin  melting 
layer  at  the  bottom  and  a  thicker  layer  at  the  top  of  the  sphere. 
They  found  that  a  significant  amount  of  melting  took  place  at  the 
upper  surface  due  to  the  significance  of  natural  convection, 
whereas  in  previous  research  [124,125],  the  effect  of  natural 
convection  on  the  melting  process  was  neglected. 

4.2.  Microencapsulated  PCM 

Micro-capsulation  is  a  kind  of  tiny  capsules  (around  0.1- 
1000  pm  in  diameter)  [127],  Micro-encapsulation  has  a  higher 
heat  transfer  rates  as  compared  to  that  of  macro-encapsulation 
[128,129]  due  to  a  substantially  higher  surface  area  to  volume 
ratio.  Meanwhile,  microencapsulation  can  tolerate  the  change  in 
volume  during  phase  change  process  [130],  Microencapsulation 
hence  is  a  prominent  technology  for  use  of  PCM  in  building 
materials.  With  the  advent  of  PCM  implemented  in  gypsum  board, 
plaster,  concrete  and/or  other  wall  covering  materials,  thermal 
storage  can  be  part  of  the  building  structure  even  for  light  weight 
buildings. 

4.2.1.  Incorporation  with  building  components 

PCMs  can  be  incorporated  with  gypsum  board,  plaster,  concrete 
and/or  other  wall  covering  materials  for  the  temperature  control 
application  purpose. 

In  general,  storage  density  and  phase  change  temperatures  are 
very  important  parameters,  since  they  determine  the  storage 
system  size  and  human  thermal  comfort.  Further,  accurate  estima¬ 
tion  of  the  PCM  enthalpy  variation  in  the  working  temperature 
range  is  essential  for  correct  mathematical  modeling  of  the  storage 
system  [131]. 

The  first  research  on  the  application  of  PCMs  in  buildings  has 
been  carried  out  by  Telkes  in  1975  [132],  Na2S04  ■  10  H20  with  a 
melting  temperature  of  32  C  was  incorporated  with  walls,  parti¬ 
tions,  ceilings  and  floors  to  serve  as  temperature  regulators. 

Determination  criteria  for  the  optimal  phase  change  temperature 
and  storage  density  were  provided  by  Zhou  [133]  as  following: 

Tm,opt  —  Tr  4-Q/htstor  (26) 

Tr  =  ttf  Td  +  tnTn/td  +  tn  (27) 

D0pt  =  tnh(Tm0pt~Tn)/pL  (28) 

where  Tmopt  is  optimal  phase  change  temperature  and  Dopt  is 
optimal  thickness  of  the  wall;  Q  represents  the  heat  absorbed  by 
unit  area  of  the  room  surface;  Tr  is  the  average  room  temperature: 
tstor  is  the  heat  storage  time;  the  subscripts  n  and  d  represent  night 
and  daytime  respectively. 

Kuznik  et  al.  [134]  developed  a  one-dimensional  model  to 
optimize  the  thickness  of  phase  change  material  only  considering 
pure  conduction.  The  simulation  was  based  on  a  lightweight  wall 
and  interior/exterior  temperature  evolutions  within  a  period  of 
24  h.  The  results  showed  that  an  optimal  PCM  thickness  exists. 
With  the  optimal  PCM  thickness  of  1  cm,  the  thermal  inertia  of  the 
building  was  doubled  and  the  PCM  can  reduce  the  temperature 
fluctuations  inside  rooms.  To  provide  guidelines  useful  in  selecting 
an  optimal  PCM  for  building  wallboard,  a  numerical  model  was 
developed  to  determine  the  optimal  melting  temperature  of  PCM 
by  Neeper  [135].  The  results  indicated  that  the  optimal  melting 
temperature  depends  on  the  average  room  temperature.  When  the 
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melting  temperature  of  PCM  closed  to  the  average  room  tempera¬ 
ture,  the  maximum  diurnal  energy  storage  was  achieved.  Xiao 
et  al.  [136]  established  a  simplified  theoretical  model  to  optimize 
an  interior  PCM  for  energy  storage  in  a  lightweight  room.  In  this 
paper  energy  conservation  equations  were  presented  to  calculate 
the  optimal  phase  change  temperature  and  the  total  amount  of 
latent  heat  capacity.  The  analytical  results  agreed  with  Neeper's 
finding  that  the  optimal  melting  temperature  depends  on  the 
average  indoor  air  temperature.  Further,  and  it  also  depends  on 
the  radiation  absorbed  by  the  PCM  panels. 

A  1-D  enthalpy  model  for  analyzing  the  thermal  performance 
of  a  shape-stabilized  PCM  floor  was  developed  by  Xu  et  al.  [137], 
By  using  the  model,  the  influence  of  various  parameters  (melting 
temperature  of  PCM,  thickness  of  PCM  layer,  latent  heat  of  fusion, 
and  thermal  conductivity  of  PCM)  on  the  thermal  performance 
was  analyzed  by  using  a  fully  implicit  the  finite-difference  method. 
It  was  found  that  for  a  given  position  or  weather  condition,  the 
suitable  melting  temperature  of  PCM  is  roughly  equal  to  the 
average  indoor  air  temperature  of  sunny  winter  days,  the  heat  of 
fusion  and  the  thermal  conductivity  of  PCM  should  be  larger  than 
120  kj /kg  and  0.5  W/(m  I<),  respectively,  and  the  thickness  of  shape 
stabilized  PCM  plate  used  under  the  floor  should  not  be  larger 
than  20  mm.  Zhou  et  al.  [138]  extended  the  investigation  of  the 
influence  of  various  parameters  from  floor  to  walls  and  ceiling  by 
using  a  verified  enthalpy  model.  The  model  was  validated  by 
experimental  work,  and  simulated  results  agreed  well  with  the 
measured  data.  The  numerical  results  showed  that  for  the  present 
conditions,  the  optimal  melting  temperature  is  about  20  C  and 
the  heat  of  fusion  should  not  be  less  than  90  kj/kg,  thin  PCM  plates 
with  large  areas  are  advantageous  and  the  effect  of  PCM  plates 
located  at  the  inner  surface  of  interior  wall  was  superior  to  that  of 
exterior  wall  (the  south  wall).  A  similar  result  to  Xu  et  al.  [137]  was 
also  obtained  for  the  condition  of  the  PCM  wall  and  ceiling  that  the 
thermal  conductivity  should  be  larger  than  0.5  W/(mK). 

Athienitis  et  al.  [139]  experimentally  and  numerically  studied 
the  application  of  PCM  in  building  envelope  components  for 
thermal  storage.  A  1-D  non-linear  enthalpy  model  was  developed 
to  simulate  the  transient  heat  transfer  process  in  the  walls  using 
explicit  finite  difference  scheme.  The  numerical  result  showed  that 
the  utilization  of  the  PCM  gypsum  board  may  reduce  the  max¬ 
imum  room  temperature  by  about  4  °C  during  the  day  time  and 
can  reduce  the  heating  load  at  night  significantly. 

The  numerical  simulation  of  passive  heating  with  PCMs  was 
carried  out  by  Chen  et  al.  [140]  by  using  a  1-D  nonlinear 
mathematical  model.  The  effective  heat  capacity  method  was  used 
to  simulate  the  PCM  heat  transfer  problem  by  using  the  software 
MATLAB.  The  result  indicated  that  applying  proper  PCM  to  the 
inner  surface  of  the  north  wall  in  the  ordinary  room  can  not  only 
enhance  the  indoor  thermal  comfort  dramatically,  but  also 
increase  the  utilization  rate  of  the  solar  radiation.  Consequently 
the  heating  energy  consuming  is  decreased  and  the  goal  of  saving 
energy  has  been  achieved.  The  energy-saving  rate  of  heating 
season  can  get  to  17%  or  higher. 


4.2.2.  Microencapsulated  phase  change  slurry  (MPCS) 

Microencapsulated  phase  change  slurry  (MPCS)  as  a  new 
technique  has  been  getting  more  and  more  attention  since  this 
technique  serves  as  both  the  heat  transfer  fluids  and  energy 
storage  media.  Consequently  it  can  increase  the  thermal  storage 
capacity  and  improve  the  energy  efficiency  of  the  thermal  system. 
Due  to  the  MPCS'  increasingly  important  role,  some  review  papers 
on  the  properties,  heat  transfer  characteristics  and  application  of 
MPCS  have  been  published  in  recent  years  [127,141-143],  As  a 
kind  of  suspension,  the  heat  transfer  and  fluid  flow  characteristics 
of  MPCS  are  different  from  those  of  the  marco-encapsulated  PCMs. 


In  accordance  to  the  review  of  Delgado  et  al.  [142],  the  heat 
transfer  characteristics  of  the  MPCS  can  be  classified  into  forced 
convection  and  natural  convection. 


4.2.2.I.  Laminar  forced  convection  heat  transfer.  The  MPCS  used  to 
enhance  convective  heat  transfer  is  generally  laminar  due  to  the  high 
viscosity.  Charunyakorn  et  al.  [144]  firstly  conducted  a  numerical 
simulation  of  MPCS  flow  in  circular  tubes  under  different  boundary 
conditions.  The  energy  equation  was  formulated  by  taking  into 
consideration  both  the  heat  absorption  (and  release)  during  phase 
transition  and  solved  using  an  implicit  finite  difference  method.  Their 
parametric  study  showed  that  the  bulk  Stefan  number  and  volumetric 
particle  concentration  are  the  two  dominant  parameters.  Their  model 
also  showed  that  the  heat  transfer  enhancement  ratio  of  the  fluid 
could  be  2-4  times  of  that  of  water.  Zhang  and  Faghri  [145]  modified 
the  model  of  Charunyakorn  [144]  to  include  the  effects  of  the 
microcapsule  walls,  initial  sub-cooling  and  the  melting  temperature 
range.  A  temperature  transforming  model  was  used  to  solve  the 
melting  in  the  microcapsule  by  instead  of  a  quasi-steady  model. 
The  numerical  results  showed  a  very  good  agreement  with  the 
experimental  results  of  Goel  et  al.  [146],  Their  numerical  results 
suggested  that  the  differences  between  the  experimental  results  of 
Goel  et  al.  and  numerical  predictions  of  Charunyakorn  could  be 
decreased  greatly  if  these  additional  factors  are  considered.  The 
authors  however  did  not  give  any  correlation  or  other  similar 
criteria  that  could  be  used  for  future  design. 

Complicated  source  terms  are  applied  to  the  theoretical  models 
above  to  account  for  the  phase  change  effects.  In  order  to  simplify 
the  simulation  process,  Alisetti  and  Roy  [147]  developed  a  simpler 
effective  specific  heat  model  to  study  the  heat  transfer  in  MPCS.  The 
results  showed  that  the  exact  form  of  the  specific  heat  function  is  not 
critical  as  long  as  the  latent  heat  is  incorporated  correctly  within  a 
finite  melting  temperature  range.  Based  on  the  concept  of  effective 
specific  heat,  Roy  and  Avanic  [148]  developed  a  model  to  analyze  the 
laminar  forced  convection  heat  transfer  to  a  PCM  suspension  in  a 
circular  duct  under  constant  wall  heat  flux  condition.  The  model  has 
been  verified  with  experimental  data.  The  results  demonstrated  that 
Stefan  number  is  the  only  parameter  that  has  a  significant  impact  on 
the  thermal  performance.  Sub-cooling  effect  only  can  affect  the  heat 
transfer  performance  at  very  low  heat  fluxes  or  when  the  inlet 
temperature  is  much  lower  than  the  melting  point.  A  simple 
correlation  to  predict  the  wall  temperature  rise  as  a  function  of  the 
tube  length  was  also  developed  for  future  design. 

Hu  and  Zhang  [149]  introduced  a  model  to  provide  a  novel  insight 
for  the  forced  convective  heat  transfer  enhancement  of  MPCS  flowing 
through  a  circular  tube  with  constant  heat  flux.  The  model  was 
developed  based  on  effective  specific  heat  capacity  method  and  used 
to  analyze  the  influence  of  various  factors  on  heat  transfer  enhance¬ 
ment.  A  modified  Nusselt  number  was  introduced  since  the  conven¬ 
tional  Nusselt  number  correlations  for  internal  flow  of  single  phase 
fluids  are  not  suitable  for  accurately  describing  the  heat  transfer 
enhancement  with  MPCS.  The  modified  local  Nusselt  number  is 
inversely  proportional  to  the  dimensionless  wall  temperature.  The 
numerical  results  indicated  that  the  exact  nature  of  the  phase  change 
process  strongly  affects  the  degree  of  convective  heat  transfer 
enhancement.  The  Stefan  number  and  the  PCM  concentration 
resulted  to  be  the  most  important  parameters  on  the  improvement 
of  heat  transfer  in  MPCS.  Zeng  et  al.  [150]  analyzed  experimentally 
and  numerically  the  enhanced  convective  heat  transfer  mechanism 
of  MPCS  in  the  thermal  fully  developed  range.  An  enthalpy  model 
was  proposed  and  developed  for  simulating  the  forced  convective 
heat  transfer  characteristics  of  laminar  flow  with  MPCS  in  a  circular 
tube  under  constant  wall  heat  flux.  The  results  showed  that  in 
the  phase  change  heat  transfer  region  the  Stefan  number  and  the 
dimensionless  phase  change  temperature  range  number  are  the 
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most  important  parameters  influencing  the  Nusselt  number  fluctua¬ 
tion  profile  and  the  dimensionless  wall  temperature.  The  bulk  fluid 
Reynolds  number,  particle  diameter  and  the  microcapsule  volumetric 
concentration  also  influence  the  Nusselt  number  profile  and  the 
dimensionless  wall  temperature  but  are  independent  of  the  phase 
change  process. 

4.22.2.  Turbulent  forced  convection  heat  transfer.  A  majority  of  past 
studies  have  considered  laminar  forced  convection  heat  transfer  to 
MPCM  and  only  a  few  numerical  investigations  on  turbulent 
forced  convection  heat  transfer  have  been  reported.  However, 
the  turbulent  flows  occur  in  majority  of  engineering  applications 
of  MPCS.  In  order  to  investigate  the  turbulent  heat  transfer  to  PCM 
suspensions,  Roy  and  Avanic  [151]  presented  an  effective  specific 
heat  capacity  model  for  PCM  suspensions  in  a  circular  tube  with 
constant  heat  flux  under  turbulent  flow  condition.  The  numerical 
results  were  found  to  agree  with  previously  published 
experimental  data  and  showed  that  this  effective  specific  heat 
model  can  effectively  model  turbulent  heat  transfer  with  PCM 
suspensions.  The  most  primary  parameter  in  heat  transfer  was  the 
Stefan  number.  The  melt  temperature  range  and  degree  of  sub¬ 
cooling  are  other  two  important  parameters.  Royon  and  Guiffant 
[152]  developed  a  numerical  model  to  describe  the  thermal 
behavior  of  a  slurry  of  PCM  slurries  flow  in  a  circular  duct  under 
different  flow  parameters  (Reynolds  number).  The  simulation 
results  demonstrated  that  the  influence  of  Reynolds  number  on 
the  minimum  length  of  the  heat  exchanger  is  relatively  weak. 

4.2.23.  Natural  convection  heat  transfer.  Inaba  et  al.  [153] 
developed  a  2-D  model  to  study  the  fluid  flow  and  heat-transfer 
characteristics  for  Rayleigh-Benard  natural  convection  of  non- 
Newtonian  PCM  slurry.  They  found  that  Rayleigh  number, 
Prandtl  number  and  aspect  ratio  could  be  the  primary  factors  for 
most  of  Newtonian  and  non-Newtonian  fluids  to  evaluate  the 
natural  convection.  A  modified  Stefan  number  was  defined  in  the 
paper  to  evaluate  the  natural  convection  in  a  PCM  slurry.  Later, 
Inaba  et  al.  [154]  presented  another  2-D  model  to  study  the 
natural  convection  heat  transfer  of  a  MPCS.  It  was  found  that  the 
natural  convection  effect  and  heat  transfer  enhancement  were  due 
to  the  contribution  of  the  latent  heat  transfer. 

4.3.  Summary  of  numerical  models 

It  can  be  seen  that  the  numerical  models  collected  above  for 
encapsulated  PCMs  all  have  their  limitations.  As  the  practical 
phase  transformations  in  different  applications  are  complicated, 
and  the  thermal  conditions  are  not  ideal.  Various  assumptions 
were  set  up  for  each  numerical  solution  according  to  different 
numerical  simulation  purposes.  There  is  none  of  numerical  models 
can  describe  the  phase  change  problem  properly.  Hence,  the 
numerical  modeling  of  LHTES  only  approximately  approaches 
but  not  accurately  presents  the  practical  thermal  performance  of 
LHTES.  Therefore,  more  investigation  should  be  carried  out  to 
study  the  phase  change  problems  so  as  to  provide  a  more  detailed 
insight  of  the  actual  heat  transfer  behaviors  inside  PCMs  during 
the  phase  transformations. 

5.  Conclusions 

This  paper  presented  a  comprehensive  review  of  mathematical 
and  numerical  methods  applied  to  the  solutions  of  phase  change 
problems.  Firstly  heat  transfer  mechanisms  during  the  charging  and 
discharging  processes  were  discussed  in  this  review.  The  heat 
transfer  mechanisms  play  the  most  important  role  on  the  numerical 
results.  Consequently,  it  is  important  to  weight  the  percentages  of 


conduction  and  convection  heat  transfer  taken  in  each  stage.  Then, 
it  presented  the  fundamental  mathematical  descriptions  of  the 
phase  change  phenomenon,  the  Stephan  problem  and  Neumann 
problem.  At  last,  a  description  of  the  numerical  solutions  for  phase 
change  problems  based  on  considering  only  conduction  and  con¬ 
sidering  also  convection.  The  PCMs  have  been  applied  in  thermal 
energy  storage  applications  widely  due  to  their  attractive  features. 
Various  numerical  models  for  encapsulated  PCMs  in  terms  of 
macro-encapsulated  and  micro-encapsulated  PCM  were  also  ana¬ 
lyzed  in  this  paper.  Nevertheless,  the  incorporation  of  PCMs  in  a 
particular  application  calls  for  an  analysis  that  will  enable  the 
researcher  to  understand  the  heat  transfer  principle  inside  PCM 
correctly.  From  this  survey,  some  further  works  are  recommended: 

•  Further  investigation  should  be  carried  out  to  obtain  adequate 
knowledge  of  the  phase  change  problems  so  as  to  minimize 
the  assumptions  made  in  the  numerical  models. 

•  Development  of  new  numerical  model  to  broaden  its  applica¬ 
tion  in  phase  change  problems  and  improve  the  accuracy  of  its 
results. 

•  Assessment  of  the  stability  and  convergence  on  the  numerical 
results  is  always  required,  and  the  numerical  predictions  should 
always  be  validated  by  using  appropriate  experimental  data. 

•  Unified  dimensionless  parameters  needs  to  be  developed  to 
make  the  gained  knowledge  applicable  to  other  cases. 

•  Life  cycle  analysis,  both  economical  and  energetic  feasibility  of 
PCM  application  should  be  performed  as  LHTES  is  a  more 
expensive  form  of  thermal  storage,  but  few  scientific  publica¬ 
tions  covered  this  matter. 
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